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Abstract Intracellular biochemical networks fluctuate dynamically due to 
various internal and external sources of fluctuation. Dissecting the fluctuation 
into biologically relevant components is important for understanding how a 
cell controls and harnesses noise and how information is transferred over 
apparently noisy intracellular networks. While substantial theoretical and 
experimental advancement on the decomposition of fluctuation was achieved 
for feedforward networks without any loop, we still lack a theoretical basis 
that can consistently extend such advancement to feedback networks. The 
main obstacle that hampers is the circulative propagation of fluctuation by 
feedback loops. In order to define the relevant quantity for the impact of 
feedback loops for fluctuation, disentanglement of the causally interlocked 
influence between the components is required. In addition, we also lack an 
approach that enables us to infer non-perturbatively the influence of the 
feedback to fluctuation as the dual reporter system does in the feedforward 
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network. In this work, we resolve these problems by extending the work on 
the fluctuation decomposition and the dual reporter system. For a single-loop 
feedback network with two components, we define feedback loop gain as the 
feedback efficiency that is consistent with the fluctuation decomposition for 
feedforward networks. Then, we clarify the relation of the feedback efficiency 
with the fluctuation propagation in an open-looped FF network. Finally, 
by extending the dual reporter system, we propose a conjugate feedback 
and feedforward system for estimating the feedback efficiency only from the 
statistics of the system non-perturbatively. 

Keywords Fluctuation • Linear Noise Approximation • Noise Decomposi¬ 
tion • Information Flow ■ Dual Reporter System 


1 Introduction 


1.1 Dissecting fluctuation in biochemical networks 


Biochemical molecules in a cell fluctuate dynamically because of the stochas¬ 
tic nature of intracellular reactions, fluctuation of the environment, and the 


spontaneous dynamics of intracellular networks 18 32 40 . Some part of the 


fluctuation is noise that impairs or disturbs the robust operation of the intra¬ 
cellular networks. The other part, however, conveys information on complex 


dynamics of various factors inside and outside of the cell 20,31 . Dissect 


ing fluctuation into distinct components with different biological roles and 
meanings is crucial for understanding the mechanisms how a cell controls 
and harnesses the noise and how information is transferred over apparently 
noisy intracellular network [5j|9]. 

The intracellular fluctuation is generated from the intracellular reactions 
or comes from the environment, and then propagates within the intracellular 
networks. The fluctuation of individual molecular species within the networks 
is therefore a consequence of the propagated fluctuation from the different 
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sources. Decomposition of the fluctuation into the contributions from the 
sources is the indispensable step for understanding their biological roles and 
relevance. When fluctuation propagates from one component to another uni- 
directionally without circulation, the fluctuation of the downstream can be 
decomposed into two contributions. One is the intrinsic part that originates 
within the pathway between the components. The other is the extrinsic part 
that propagates from the upstream component. Such decomposition can eas¬ 
ily be extended for the network with cascading or branching structures in 
which no feedback exists. This fact drove the intensive anatomical analysis 
of the intracellular fluctuation in the last decade. 


1.2 Decomposition of sources of fluctuation 

In order to dissect fluctuation into different components, two major strategies 
have been developed. One is to use the dependency of each component on 
different kinetic parameters in the network. By employing theoretical predic¬ 
tions on such dependency, we can estimate the relative contributions of dif¬ 
ferent components from single-cell experiments with perturbations. Possible 
decompositions of the fluctuation were investigated theoretically for various 
networks such as single gene expression |3 26, 27], signal transduction path¬ 
ways (4l]47 , and cascading reactions |45| . Some of them were experimentally 
tested [3 26,28 . 


The other is the dual reporter system in which we simultaneously measure 
a target molecule with its replica obtained by synthetically duplicating the 
target. From the statistics of the target and the replica, i.e, mean, variance, 
and covariance, we can discriminate the intrinsic and extrinsic contributions 
to the fluctuation because the former is independent between the target and 
the replica whereas the latter is common to them. The idea of this strat¬ 


egy was proposed and developed in 27 43 , and verified experimentally for 
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different species 11 23 33 . Its applicability and generality were further ex¬ 


tended 4| 6 16 17 35 


Now these strategies play the fundamental role to design single-cell ex¬ 
periments and to derive information on the anatomy of fluctuation from the 
experimental observations [8[|l5}[29,34p7[|44]. 

1.3 Feedback regulation and its efficiency 


Even with the theoretical and the experimental advancement in decomposing 
fluctuation, most of works focused on the feedforward (FF) networks in which 
no feedback and circulation exist. As commonly known in the control theory 
[7 39 , feedback (FB) loops substantially affect fluctuation of a network by 
either suppressing or amplifying it. Actually, the suppression of fluctuation 
in a single gene expression with a FB loop was experimentally tested in [2j 
earlier than the decomposition of fluctuation. While the qualitative and the 
quantitative impacts of the FB loops were investigated both theoretically 
and experimentally [I| [2H[24{[25l[35|[42|[46] , we still lack a theoretical basis 
that can consistently integrate such knowledge with that on the fluctuation 
decomposition developed for the FF networks. 

The main problem that hampers the integration is the circulation of fluc¬ 
tuation in the FB network. Because fluctuation generated at a molecular 
component propagates the network back to itself, we need to disentangle 
the causally interlocked influence between the components to define the rel¬ 
evant quantity for the impact of the FB loops. From the experimental point 
of view, in addition, quantification of the impact of FB loops by perturba¬ 
tive experiments is not perfectly reliable because artificial blocking of the 
FB loops inevitably accompanies the change not only in fluctuation but also 
in the average level of the molecular components involved in the loops. It is 
quite demanding and almost impossible for most cases to inhibit the loops by 
keeping the average level unchanged. We still lack an approach that enables 
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us to infer the influence of the FB non-perturbatively as the dual reporter 
system does. 


1.4 Outline of this work 


In this work, we resolve these problems by extending the work on the fluctu¬ 


ation decomposition 27 41 and the dual reporter system 11 43 . By using 


a single-loop FB network with two components and its linear noise approxi¬ 
mation (LNA) |10, T9p6] , we first provide a definition of the FB loop gain as 
FB efficiency that is consistent with the fluctuation decomposition in [27,41 


Then, we clarify the relation of the FB efficiency with the fluctuation prop¬ 
agation in a corresponding open-looped FF network. Finally, by extending 
the dual reporter system, we propose a conjugate FB and FF system for 
estimating the feedback efficiency only from the statistics of the system non- 
perturbatively. We also give a fluctuation relation among the statistics that 
may be used to check the validity of the LNA for a given network. 

The rest of this paper is organized as follows. In Sect. [2j we review the 
decomposition of fluctuation for a simple FF system derived in 27j[4l] by 
using the LNA. In Sect. [3j we extend the result shown in Sect. [2] to a FB 
network by deriving a decomposition of the fluctuation with feedback. Using 
this decomposition, we define the FB loop gain that is relevant for quantify¬ 
ing the impact of the FB to the fluctuation. In Sect. [4j we give a quantitative 
relation of the loop gain in the FB network with the fluctuation propagation 
in a corresponding open-looped FF network. In Sect. [5j we propose a conju¬ 
gate FF and FB network as a natural extension of the dual reporter systems 
used mainly for the FF networks. We clarify that the loop gain can be es¬ 
timated only from the statistics, i.e., mean, variances, and covariances, of 
the conjugate network. We also show that a fluctuation relation holds among 
the statistics, which generalizes the relation used in the dual reporter sys- 
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tem. In Sect. [6j we discuss a link of the conjugate network with the directed 
information, and give future directions of our work. 


2 Fluctuation Decomposition and Propagation in a Small 
Biochemical Network 


In this section, we summarize the result for the decomposition of fluctuation 


obtained in 27 41 by using the LNA, and also its relation with the dual 


reporter system that was employed in ll,|23|33j to quantify the intrinsic 
and extrinsic contributions from the experimental measurements. 


2.1 Stochastic chemical reaction and its linear noise approximation 

Let us consider a chemical reaction network consisting of N different molecu¬ 
lar species and M different reactions. We assume that the stochastic dynamics 
of the network is modeled by the following chemical master equation: 


n ~ s k) - a k (n)F(t, n )], 


(1) 


k =1 


where n = (m,..., unj E N> 0 is the numbers of the molecular species, 
P(t, n) is the probability that the number of molecular species is n at t, 
and Ofc(n) eM55 0 and are the propensity function and the stoichiometric 


vector of the fcth reaction, respectively 12 14 19 . The propensity function 


characterizes the probability of occurrence of the fcth reaction when the num¬ 
ber of the molecular species is n, and the stoichiometric vector defines the 
change in the number of the molecular species when the fcth reaction occurs. 

In general, it is almost impossible to directly solve Eq. both analyti¬ 
cally and numerically because it is a high-dimensional or infinite-dimensional 
differential equation. To obtain insights for the dynamics of the reaction net¬ 


work, several approximations have been introduced 12 19 . Among others, 
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the first-order approximation is the deterministic reaction equation that is 
described for the given propensity functions and the stoichiometric vector as 


df2u(t) 
d t 


= Sa(f2u(t )), 


( 2 ) 


where h? is the system size, u £ M> 0 is the concentration of n as u = n/£2, 
a(n) := (ai(n)),..., aM(n)) T , and S+ t k := Sk is the stoichiometric matrix. 
Equation § was successfully applied for chemical reaction networks with 
large system size where the fluctuation of the concentration of the molecu¬ 
lar species can be neglected. When the system size is not sufficiently large, 
however, Eq. ([2]) is not relevant for analyzing the fluctuation of the network. 
The LNA is a kind of the second-order approximation of the Eq. 0 that 
characterizes the fluctuation around a fixed point, n, of Eq. ([2]) that satis¬ 
fies Sa(n) = 0. The stationary fluctuation of the network around n is then 
obtained by solving the following Lyapunov equation 10,19,46 


(K(n)E) + ( K(n)E) r + D{n ) = 0, 


( 3 ) 


where 


v , n_ dSa(n) _. / \ 

-***yTl J .— Qji ■ ’ j Sj,kSj,kQ j k\'H') ? 


( 4 ) 


and E is the covariance matrix of n. When the propensity function a(n) 
is affine with respect to n, the dynamics of Qu(t ) determined by Eq. 0 is 
identical to that of the first cumulant of n, i.e., the average of n, as 


d (™) 

d t 


= Sa((n)), 


( 5 ) 


where (n) := n¥(t, n). In addition, the second cumulant, i.e., the covari¬ 

ance matrix, follows the Lyapunov equation as 


d E 
d t 


= (K((n))E) + (. K((n))Ef + £>((«)). 


( 6 ) 










Therefore, if the propensity function a(n) is affine, the stationary fluctuation 
of n is exactly described by Eq. ([3]). For a non-affine a(n), Eq. ([ 5 ]) and Eq. ([h]) 
can also be regarded as an approximation of the full cumulant equations 


by the cumulant closure 12 under which we ignore the influence of the 


second and the higher order cumulants to Eq. ([5]), and that of the third 
and the higher order ones to Eq. (JbJQ Even though the propensity function 
a(n) is not affine, Eq. ^ (or Eq. §) can produce a good approximation 
of the fluctuation, provided that the fixed point, n, is a good approximation 
of the average, (ra), and that the local dynamics around the fixed point is 
approximated enough by its linearization. In addition, compared with other 
approximations, the LNA enables us to obtain an analytic representation of 
the fluctuation because Eq. Q is a linear algebraic equation with respect 
to S. Owing to this property, the LNA and its variations played the crucial 
role to reveal the analytic representation of the fluctuation decomposition in 
biochemical reaction networks |27j46| . As in these previous works, we employ 
the LNA to obtain an analytic representation for the feedback efficiency. 


2.2 Decomposition of fluctuation 

Starting from the LNA, Paulsson derived an analytic result on how noise is 
determined in a FF network with two components (Fig. a (A)) [27], Here, 


we briefly summarize the result derived in 27 . Let n\ = x and ri 2 = y for 


notational simplicity, and consider the FF reaction network (Fig. [I] (A)) with 
the following propensity function and stoichiometric matrix, 


a(x, y) = (a+ (x), a x (x), a+ (x, y),a y ( x , y)), 


S = 


+ 1-10 0 
0 0 + 1 - 1 , 


(7) 


1 We prefer this interpretation of the Lyapunov equation because the LNA has 
been applied for various intracellular networks whose system size is not sufficiently 
large. 
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Fig. 1 (A) The structure of the t wo-c omponent FF network. Interpretations of this 
network as single-gene expression 11122||26|27 , two-gene regulation [46], and signal 
transduction [47| are shown. (B) The structure of the dual reporter system [ill 
(C) A schematic diagram of the FF network for two-gene regulation.(U) 
A schematic diagram of the dual reporter network for the two-gene regulation. 


Because a x depends only on x, x regulates y unidirectionally. Then, for a 
fixed point (x,y) of Eq. |2]), that satisfies 


a+{x)=a x (x), a y (x, y) = a y (x, y), 


( 8 ) 


K and D in Eq. © becomes 


K = 



H x x l T~X 

~I Hyx/ T y —Hyy/'y 



(2a x 0 \ 

ii 

V 0 2 dy 


( 9 ) 


where a x := a x (x) = a x (x), a v = a^{x 1 y) = a y (x : y), and Hi j := 


In a i I ay 


3 In j 


(x,v) 


for i,j £ {x,y}. d x and d y are the minus of the diagonal terms of K and 
represent the effective degradation rates of x and y. k yx is the off-diagonal 
term of K that represents the interaction from x to y. 


Hij is the susceptibility of the *th component to the perturbation of 
the jth one. t x := x/a x and r y := y/d v are the effective life-time of x 
and y, respectively. Except this section, we mainly use ds and ks as the 
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representation of the parameters rather than Hs and rs introduced in 27[ 
By solving Eq. ([3]) analytically, the following fluctuation-dissipation relation 
was derived in [27 as 



(Hi) 



rf • (10) 


This representation measures the intensity of the fluctuation by the coefficient 
of variation (CV)[^] and describes how the fluctuation generates and propa¬ 
gates within the network. (I) is the intrinsic fluctuation of x that originates 
from the stochastic birth and death of a;. 1/5 reflects the Poissonian nature 
of the stochastic birth and death, and 1/H XX is the effect of auto-regulatory 
FB. Similarly, (II) is the intrinsic fluctuation of y that originates from the 
stochastic birth and death of y. (Ill), on the other hand, corresponds to 
the extrinsic contributions to the fluctuation of y due to the fluctuation of 
x. The term (III) is further decomposed into (i), (ii), and (in), (i) is the 
fluctuation of x , and therefore, identical to (I), (ii) and (Hi) determine the 
efficiency of the propagation of the fluctuation from x to y, which correspond 
to the time-averaging and sensitivity of the pathway from x to y , respectively. 
This representation captures the important difference of the intrinsic and the 
extrinsic fluctation such that the intrinsic one, the term (II), can be reduced 
by increasing the average of y whereas the extrinsic one, the term (III), 
cannot. 


While Eq. (101 provides an useful interpretation on how the fluctuation 
propagates in the FF network, it is not appropriate for the extension to the 
FB network because the contribution of x to the fluctuation of y is described 
by the CV of x as the term (i). Because the fluctuation of x and y depend 


3 The former gives a notationally simpler result. 

3 The CV is defined by the ratio of the standard deviation to the mean as a x / {x }. 
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mutually if we have a FB between x and y, we need to characterize the 
fluctuation of y without directly using the fluctuation of x. To this end, we 
adopt the variances and covariances as the measure of the fluctuation and 
use the following decomposition of the fluctuation for the FF network: 


(«) (i) 


at = G x 


<T y ~ Gyy 


(I) 


C«) (i) (***) (**) (i) 

V + Gy X G x ,x x 
(//) (///) 


( 11 ) 


where we define G xx , G vv , and G yx as 


G xx — 


H :r 


Gyy ~ „ 
n, 


yy 


h 2 

G = yx 

VX (d x + dy ) dy ' 


( 12 ) 


The terms (I), (//), and (///) in Eq. (11) correspond to those in Eq. (10). 
The interpretation of the terms within (/), (II), and (III) is, however, dif¬ 


ferent. In Eq. (11), the terms (I) are interpreted as the fluctuation purely 
generated by the birth and death reactions of x and y by neglecting any con¬ 
tribution of the auto-FBs. Because the simple birth and death of a molecular 
species without any regulation follow the Poissonian statistics, the intensity 
of the fluctuation is equal to the means of the species. Thus, x « (x) and 
y ~ (y) in the terms (I) represent the generation of the fluctuation by the 
birth and death of x and y, respectively. The fluctuation generated is then 
amplified or suppressed by the auto-regulatory FBs. G xx and G vy in the 
terms (ii) account for this influence, and are denoted as auto-FB gains in 
this work. Finally, G yx in the term (in) quantifies the efficiency of the prop¬ 
agation of the fluctuation from x to y. We denote G yx as the path gain from 


x to y. If we use the notation in Eq. (10), G yx is described as 

Hyy/ T y 


Gyx — 


yHyx/Ty 

^Hyy / Ty J Hyy/Ty T H XX /T X 


(13) 


The decomposition of the fluctuation of y into (II) and (III) is consis¬ 


tent between Eq. (10) and Eq. (11) while the further decompositions within 














12 


(II) and (III) are different. It was recently revealed in [4.6,16,17 that the 
decomposition into (II) and (III) is linked to the variance decomposition 
formula in statistics as 


<4 = V\E\y(t)\X(t)]] +E[V[y(t)|A’(t)]], 


(14) 


(ii) 


(in) 


where X(t) := {x(t)-,t £ [0,f]} is the history of x(t), and E and V are the 
expectation and the variance, respectivelj0 Because the variance decompo¬ 


sition formula holds generally, Eq. (|14[) is more fundamental than Eq. (10) 


and Eq. (11) as the decomposition of fluctuation. This decomposition was 


further analyzed in 4[ 6 16 


2.3 Dual reporter system 


The decomposition, Eq. (10) or Eq. © , guides us how to evaluate the intrin¬ 
sic and the extrinsic fluctuation in y experimentally. When we can externally 
control the mean of y without affecting the term (III), we can estimate 
the relative contributions of (II) and (III) as the intrinsic and the extrin¬ 
sic fluctuation by plotting <jy as a function of the mean of y. When x and 
y correspond to a mRNA and a protein in the single gene expression, the 


translation rate works as such a control parameter 22 . This approach was 


intensively employed to estimate the efficiency of the fluctuation propagation 
in various intracellular networks 3 26 28 . 


Another way to quantify the intrinsic and the extrinsic fluctuation is the 
dual reporter system adopted in [ITp3p3l whose network structure is shown 
in Fig. 0(B0 In the dual reporter system, a replica of y is attached to the 
downstream of x as in Fig. [l] (B) where y' denotes the molecular species 


The correspondence of Eq. (1141) with Eq. (10) or Eq. (Ill is valid only when 


ay is decomposed under the conditioning with respect to the history of x(t), X(t), 
rather than the instantaneous state of x(t) at t. 

5 The dual reporter system is also called a conjugate reporter system 4,6. 
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of the replica. The replica, y ', must have the same kinetics as y, and must 
be measured simultaneously with y. If y is a protein whose expression is 
regulated by another protein, x, as in Fig. |T] (C), then y' can be synthetically 
constructed by duplicating the gene of y and attaching fluorescent probes 
with different colors to y and y' as in Fig. a (D) [II]. Under the LNA, the 
covariance between y and y' can be described as 


a 


(15) 


Thus, by using only the statistics of the dual reporter system, the intrinsic 
and the extrinsic components in y can be estimated as 


a y ~ Gyy 


(y) 


'y,y * 


(16) 


(ii) 


(in) 


where it is unnecessary to control any kinetic parameters externally. This 
approach is generalized as the conjugate reporter system together with the 


more general decomposition, Eq. (14), in 4.6 16 


3 Feedback Loop Gain in a Small Biochemical Network 


While we had substantial progress in the decomposition of the fluctuation 
and its experimental measurement for the FF networks in the last decade, 
its extension to FB networks is yet to be achieved. In this section, we extend 
the decomposition of the fluctuation in the simple FF network (Eq. ©) to 
the corresponding FB one depicted in Fig. [2] (A) . As in Eq. ([9]), I\ and D 
in the Lyapunov equation (Eq. @) for the FB network can be described as 


K = 




^ H X y/T X 

~Hyy/ T y 



( 17 ) 
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Fig. 2 (A) The structure of the two-component FB network. (B) A schematic 
diagram of the FB network for two-gene regulation. 


By defining the path gain from y back to x as 


fi x v 

xy ' d x {d x + d y y 


(18) 


we can derive a decomposition of the fluctuation of x and y as 


(I) 

rr 2 — 1 ~ Ly _i_ 

a x — i T T ^xx x T 

1 ±J X ±Jy 

2 _ 1 ^x 

ay ~ i - L r - l„ 


1 L x Ly 


(IV) 

G X yG yy y, 


GyyV,- 

(II) 


. J . Gy X G XX X, 

1 L/ x L/y \ - v - y 

(III) 


(19) 


where we define 


L x \ — 


k k 

rX/ xy rv yx 

dx (d x + dy 


LV (d, 


k k 

rx, xy rx 'y x 


dy)dy 


( 20 ) 


If the FB from y to x does not exist, i.e., k xy = 0, then L x = L y = G xy = 0 


and Eq. (19) is reduced to Eq. (11). Thus, L x and L y account for the effect 


of the FB. L x and L y are denoted as the FB loop gains in this work. Eq. (19) 


clearly demonstrates that the representation in Eq. (101 does not work with 


the FB because a 2 cannot be described with a 2 any longer. In contrast, we 
can interpret the terms in Eq. (19) consistently with those in Eq. because 


all the terms, (J), (II) and (HI), are unchanged in Eq. (19). The new term, 


(IV), in the expression for a 2 appears to account for the propagation of the 
fluctuation generated by the birth and death events of y back to x. 
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Equations (19) and (20) indicate how the FB affects the fluctuation of 
x and y. First, the fluctuation is suppressed when L x and L y are negative 
whereas it is amplified when they are positive. When d x and d y are positiv^J 
the sign of L x and L y are determined by the sign of k xy k yx . Thus, the FB 
loop is negative when x regulates y positively and y does x negatively or vise 
versa. This is consistent with the normal definition of the sign of a FB loop. 
Second, the efficiency of the FB depends on the source of the fluctuation. 
For example, when L x <C L y , e.g., the time-scale of x is much faster than y 1 


then Eq. (19) can be approximated as 


(i) 


(IV) 


r'i = G, 


1 - L„ 


1 -L„ 


GyyV + G yxG xx ^X 

■' S '—V-' 

L (II) (III) J 


Thus, the FB does not work for the term (J) that is the part of fluctuation of 
x whose origin is the birth and death events of x itself. This result reflects the 
fact that the slow FB from x to itself via y cannot affect the fast component 
of the fluctuation of aQ Finally, when L x and L y satisfy 1 — L x — L y = 
1 — k xy k yx /d x d y = 0, the fluctuation of both x and y diverges due to the 
FB. Since this condition means that the determinant of K becomes 0 and K 
is the Jacobian matrix of Eq. §> the fluctuation of x and y diverges due to 
the destabilization of the fixed point, fi, by the FB. 


6 For biologically relevant situations, d x and d y are positive because they can be 
regarded as the effective degradation rates. 

7 Note that the term (III) is affected by the FB even though its origin is the 
fast birth and death events of x. This can be explained as follows. In general, 
the path gain from x to y, G yx , becomes very small compared to others when x 
has much faster time scale than y. Thus, the term (III) becomes quite small and 
little fluctuation propagates from x to y because of the averaging effect of the slow 
dynamics of y. (Ill) represents, therefore, the slow component in the fluctuation 
of the birth and death of x that has the comparative timescale as that of y. This 
is why the slow FB can affect the term (III)- 
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No interaction 


Fig. 3 (A) The structure of the the opened FF network obtained by replicating y in 
the FB network. (B) A schematic diagram of the network in (A) for gene regulation. 
(C) The structure of the the opened FF network obtained by replicating x in the 
FB network. (D) A schematic diagram of the network in (C) for gene regulation. 


4 Relation between Fluctuation Propagation and Feedback Gain 


As shown in Sect. [3j L x and L y are quantitatively related to the efficiency of 
the FB. Because L x L y = G xy G yx holds, the loop gains are also linked to the 
propagation of the fluctuation from x to y and from y back to x. However, 
the meaning of the individual gains, L x and L y , is still ambiguous. 

To clarify the relation between the loop gains and the propagation of 
fluctuation, we consider a three-component FF network shown in Fig. 1(A) 
and (B) that are obtained by opening the FB network in Fig. 1(A). In the 
network shown in Fig. 1(A), x is regulated not by y but by its replica denoted 
as y'. We assume that y' is not driven by x, and thereby, y' —» x —> y forms 
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a FF network. K and D in Eq. for this network become 


1 -d y , 0 0 N 

ky' x d x 0 

\ 0 k x y d,y J 



( 21 ) 


Because y' is the replica of y, we assume that all the kinetic parameters of 
y' are equal to those of y as k y ' x = k yx , d y > = d y , and a y ' = a y . 


By solving Eq. ([3]) , we can obtain the following decomposition of the 
fluctuation of y as 


where 


(IV) 



(J — G xx X-\~G xy 'G 


xy'^Jy'y'y , 


(Ty G yy y T Gy X G XX X F \Gyy' Gy X y '] Gy'y/y , 


(II) 


(III) 


(V) 


( 22 ) 


G X y< — G 


xy, 


G 


yxy' 


• Gy X G X y' Gy X G X y L X L y , Gyy' • L y / *2. ( 23 ) 


The term (IV) for x is similar to the propagation of the fluctuation from y 
to x in the FB network, but it represents the propagation of the fluctuation 
from the replica y' in this opened FF network. The new term ( V ) accounts 
for the propagation of the fluctuation from y' down to y. The gain of this 
propagation has two terms, G vxy ' and G yy '. The first term, G yxy /, is the total 
gain of the fluctuation propagation from y' to x and from x to y' because 
G yxy ' = G yx G xy ' holds. In order to see the meaning of the second term G yy >, 


we need to rearrange Eq. (22) as 


<J y ~ Gyyy + G yx a x + Gyy'Gy'y'y . 


( 24 ) 
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This representation clarifies that G yy ’ describes the propagation of the fluc¬ 
tuation from y' to y that cannot be reflected to the fluctuation of the inter¬ 
mediate component, x. In addition, by solving Eq. ([3]), we can see that the 
gain Gyy> is directly related to the covariance between y and y' as 

<v,i/ = V G iw'y 1 = ^%y'- ( 25 ) 


This implies that we have at least two types of propagation of the fluctua¬ 
tion. One described by G yxy ' is that the fluctuation of the upstream, i.e., y ', 
is absorbed by the intermediate component, i.e., x, and then the absorbed 
fluctuation propagates to the downstream, i.e., y. This component does not 
convey the information of the upstream because that does not affect the co- 
variance between the upstream and the downstream. The other described by 
G yy i is that the fluctuation of the upstream propagates to the downstream 
without affecting the intermediate component. This fluctuation conveys the 
information on the upstream to the downstream because it is directly linked 
to their covariance. The fact that L y is related to the latter indicates that 
the FB efficiency is directly linked to the information transfer in the opened 
loop from y' to y. By considering another opened loop network where the 
replica of x is introduced as in Fig. [3] (C), we can obtain the following result: 


(ii) 


(in) 


n y — G yy y + G yx iG x ' x 'X , 

o' x G xx x -\~ G xy G yy y T \G XX > ~\- G xyx '] G x > x >x , 


(i) 


(IV) 


(vi) 


(26) 


where G xx j \= L 2 X / 2, and we also have 


® X 1 ,X 


sfG x 



(27) 


By combining Eq. (25) and Eq. (27), we can estimate the loop gains, L x and 


L y , by only measuring the averages and covariances of the opened networks 
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as follows: 


L X 

2 


&x' 


( x ') ’ 


a y',v 

W) 


(28) 


While this strategy sounds eligible theoretically, it has an experimental diffi¬ 
culty in constructing the opened networks. In the opened network, the replica, 
e.g., y ', must be designed so that it is free from the regulation of x by keeping 
all the other properties and kinetic parameters the same as those of y. For 
example, if x and y are regulatory proteins and if they regulate each other as 
transcription factors as Fig. [2] (B), the replica, y', must not be regulated by x 
but its expression rate must be equal to the average expression rate of y under 
the regulation of x as in Fig.[3](B). This requires fine-tuning of the expression 
rate of y' by modifying the DNA sequences relevant for the rate. In order to 
conduct this tuning, we have to measure several kinetic parameters of the 
original FB networks that undermines the advantage of the opened network 
that measurements of the kinetic parameters are unnecessary to estimate L x 


and L y via Eq. (28). 


5 Estimation of Feedback Loop Gain by a Conjugate FB &; FF 
Network 

As a more promising strategy for the measurement of the loop gains, we pro¬ 
pose a conjugate FB and FF network that is an extension of the dual reporter 
system for the estimation of the intrinsic and the extrinsic components. In 
the conjugate network, we couple the original FB network with a replica that 
is opened as in Fig. [4] (A), x and y are the same as the original FB network. 
The replica, y', is regulated by x as y is but does not regulate x back. Thus, 
x and the replica, y', form an FF network. If x and y are regulatory proteins 
as in Fig. [2] (B), the replica y' can be engineered by duplicating the gene 
y with its promoter site, and by modifying the coding region of the replica 
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(A) (B) 
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OO 
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No interaction 
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Fig. 4 (A) The structure of the conjugate FB and FF network obtained by repli¬ 
cating y in the FB network. (B) A schematic diagram of the network in (A) for 
gene regulation. (C) The structure of the conjugate FB and FF network obtained 
by replicating x in the FB network. (D) A schematic diagram of the network in 
(C) for gene regulation. 


so that y' looses the affinity for binding to the regulatory region of x as in 
Fig. 0 (B). This engineering is much easier than that required for designing 
the opened network. 


Next, we show how to use this conjugate network to measure the loop 
gains. For this network, K and D in Eq. © become 



—d k 0 \ 

^xy 


'ct x 0 O' 

K = 

O 

1 

, D = 

0 dy 0 


iky' X 0 dy' j 


o 

o 

©I 


(29) 


Because the replica y' affects neither x nor y , the fluctuation of x and that 


of y are the same as those of the FB network in Eq. (19). The variance of 
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the replica y' can be decomposed as 



y'x K - T xx 


(1 -L x -L y )(l-L y> ) 

(Gy'y + Gy' X y)Gy y y, 


(1 L x Ly)( 1 Ly') ■ 


(V) 


(30) 


where 


Lyl — 


/~i _ /-i 

r> i '-'y'y' ~ '- z ; 


yyi 


Gy>x Gyx) 


Gy'y — G 


yy' > 


Gy'xy - Gy X G X y • 


(31) 


Rearrangement of this equation leads to 


2 _ 2 

" 1 - Ly/ 2 


GyyV- 


(32) 


In addition, we have the following expression for the covariance between y 
and y' as 


1 


&v.v' — 


(1 L x Ly) 

_ 2_ 1 

y 1-V2 


G ra G 

XX x + 


(1 + L X )Ly 


(1 L x Ly)( 1 Ly') 


GyyV, 


(33) 


GyyV- 


A similar result can be obtained by replicating x as in Fig. [4] (C) and (D). 


By using these relations, we have 


(a) 


2 

-&X 

&x,x 

' -O'* 

a 2 y , 


a y,y‘ 




( 

» 


(b) (c) 

/-^'-*-■> 

1 _ / G xx \ 

1/2 + G XX /{F X - F X ') ~ 2 V " F x - F x ,x' ) ’ 

_I_ = 2 f l_^ 

1/2 + Gyy/(Fy ~ Fy'), \ Fy ~ Fy y' J 

S " V* ^ ^ V* ^ 

(b) (C) 


(34) 


where F x := er^/ (x), := a 2 / (y), F XtX ' ■= a x ,x>/ {x), and F y>y i := a y<y > / {y) 

are the Fano factors of x and y and normalized covariances. This result 
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Fig. 5 (A) A schematic diagram of the conjugate FB and FF network used for the 
numerical simulation. (B) The distributions oip(x,y) ,p(x,y'), and p(x',y) sampled 
from the simulation for different parameter values of K x . The other parameters are 


fo = go= 400, K y = 100, 


distributions ofp(x, 
the nullclines of Eq. 
with the relations, (aj 


and d x = d y = 1. The blue, 


red, and green dots are the 


, p(x, y'), and^a/, y), respectively. Solid and dashed lines are 
d2|) defined as da:/dt = 0 an d dy/dt = 0 . (C) ( L x , L y ) estimated 
(b), and (c), in Eq. ( |34| | from the numerical simulation. The 
parameter values used were the s ame as in (B). For each estimation, the means and 
the covariances required in Eq. (341 were calculated from 9 x 10 4 samplings. For 
each parameter value, we calculated the estimates 10 times to see their variation. 
Black markers show the analytically obtained values of (L x , L y ) for each parameter 
value. (D) A plot of the estimates of F x + F x i — 2F x , x i and F y +F y i — 2 F y y derived 
in Eq. (|35|). The parameter values used were the same as those in (C). 


indicate that we have multiple ways, (a), (b), and (c), to estimate L x and 
L y from the statistics of the conjugate network. In addition, we also have a 
fluctuation relation that holds between the statistics as 


F A- F , — 2F / — 2(7 

z X I x x' x,x' — ^ y ~ T X 


F„ 


Fyl "dFy^y! — 2Gyy, 


( 35 ) 


which is the generalization of Eq. (16) for the dual reporter system. 
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5.1 Verification of the relations by numerical simulation 


We verify Eq. (34) and Eq. (35) by using numerical simulation. For the sim¬ 


ulation, we use a conjugate network of gene regulation in which the replicas 
of both x and y are involved to measure L x and L y , simultaneously as in 
Fig. [ 5 J A). For the variable n = (x, y, x' , y') T , the stoichiometric matrix and 
the propensity function are 

( f(y ) \ 

-d x x 

9(x) 

—d v y 

f(y) 

—d x x' 

9 0) 

\~dyV 7 

The simulation is conducted by the Gillespie’s next reaction algorithm 13]. 

First, we test linear negative feedback regulation defined by f(y) := max f 0 [1 — ^-], 0 


S = 


/1-1 0 0 0 0 0 0 ^ 
0 0 1 -1 0 0 0 0 
0 0 0 0 1 -1 0 0 
\0 0 0 0 0 0 1 - 1 / 


l(x) = 


(36) 


and g(x) := max 


9o±, 0 


under which the LNA holds exactly as long as its 


trajectory has sufficiently small probability to reach the boundaries of x = 0 
or y = 0. In Fig. [5](B) , the distributions of p(x,y), p(x,y'), and p(x',y) are 
plotted for different parameter values of K x . The FB is strong for small K x 


whereas it is weak for large K x . (L x . L y ) estimated by Eq. (341 for the param¬ 


eter values are plotted in Fig. [5](C). The three estimators, (a), (b), and (c) in 


Eq. (34), are used for comparison. For this simulation, all the estimators work 


well, but they have slightly larger variability in L x compared with that in 
L y for large values of |La;| and \L y \. In addition, when \L X \ and \L y \ are very 
small, i.e. much less than 1, the estimators show relatively larger variability 
and bias, suggesting that the estimation of very weak FB efficiency requires 
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Fig. 6 (A) The distributions of p(x,y) , p(x,y'), and p(x',y) sampled from the 

simulation for different parameter values of n x and n y where n x = n y = |n|. 
The other parameters are fo = go = 300, K x = K y = 150, and d x = d y = 1. 
The blue, red, and green dots are the distributions of p(x, y) .^p(x, y'), and p(x', y), 
respectively. Solid and dashed curves are the nullclines of Eq. J2J) defined as dx/dt = 
0 and dy/dt = 0 . (B) ( L x ,L y ) estimated with the relations, (a), (b), and (c), in 
Eq. ( |34[ ) from the numerical simulation. The same parameter values were used as 
in (Aj. For each estimation, the means and the covariances required in Eq. (341 
were calculated from 9 x 10 4 samplings. For each parameter value, we calculated 
the estimates 10 times to see their variation. Black markers show the analytically 
obtained values of ( L x ,L y ) for each parameter value. ( C) A plot of the estimates 
of F x + F x / — 2 F x y and F y + F y i — 2 F y y i shown in Eq. (351. The parameter values 
used were the same as those in (B). 


more sampling. For the same parameter values, we also check Eq. (351 in 
Fig.[5}D). Both F x + F x i — 2 F XtX > and F y + F y r — 2F y y localize near 2 ir¬ 
respective of the parameter values. As Fig. [5jC) and (D) demonstrate, the 
estimators obtained from the simulations agree with the analytical values of 
(L x ,L y ), and the fluctuation relation also holds very robustly. 


In order to test how nonlinearity affects the estimation, we also investi¬ 
gated non-linear negative feedback regulation defined by f(y) := /o ) n y ; 

and g(x) := ffo i +/J/k )"* • change the Hill coefficients, n x and n y , by 
keeping the fixed point unchanged as in Fig. [djjA). Compared with the linear 
case, the variability of the estimators is almost similar even though the feed¬ 
back regulation is nonlinear (Fig. |6|B)). In addition, the estimators show a 
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good agreement with the analytical value except for very large value of \L X 


and \L y \. This suggests that Eq. (34) works as good estimators when the 
trajectories of the system are localized sufficiently near the fixed point as in 
Fig.gA). For the very large value of \L X \ and \L y \ where |n| = 2 7 / 2 , all the 


estimators are slightly biased towards smaller values, and (b) in Eq. (34) has 
larger variance than the others. In addition, similarly to the linear case, the 
estimators require larger sampling when \L X \ and \L y \ are much less than 1. 
Even with the nonlinear regulation, the fluctuation relation holds robustly 
as shown in Fig. §C), which is consistent with the good agreement of the 
estimators of L x and L y with their analytical values as in Fig. [6^B) . 

All these results indicate that the estimators obtained in Eq. ( |34| ) can be 
used to estimate the FB efficiency as long as the efficiency is moderate and 
the trajectories of the system are localized near the fixed point. 


6 Discussion 

In this work, we extended the fluctuation decomposition formula obtained 
for the FF network to the FB network. In this extension, the FB loop gains 
are naturally derived as the measure to quantify the efficiency of the FB. By 
considering the opened FF network obtained by opening the loop of the FB 
network, the relation between the loop gains and the fluctuation propagation 
in the FF network was clarified. In addition, we proposed the conjugate FB 
and FF network as a methodology to quantify the loop gains by showing 
that the loop gains are estimated only from the statistics of the conjugate 
network. By using numerical simulation, we demonstrated that the loop gains 
can actually be estimated by the conjugate network while we need more 
investigation on the bias and variance of the estimators. Furthermore, the 
fluctuation relation that holds in the conjugate network was also verified. 
We think that our work gives a theoretical basis for the conjugate network 
as a scheme for experimental estimation of the FB loop gains. 





26 


As for the further problems on the efficiency of FB and its quantification, 
the generality of the FB loop gain should be clarified. In the case of the 
FF network, the fluctuation decomposition proposed initially by using the 
LNA was generalized as the variance decomposition formula with respect 
to the conditioning of the history of the upstream fluctuation [4,17 . In the 
case of the FB network, similarly, the decomposition and the FB loop gains 
were obtained and defined via the LNA. However, its generalization is not 
straightforward because the contribution of the interlocked components can¬ 
not be dissected in the FB network because of the circulative flow of the 
fluctuation. The previous work on the relation between the fluctuation de¬ 
composition and the information transfer for the FF network [1] suggests 
that the FB loop gains can also be interpreted in terms of the information 
transfer. As a candidate for such information measure, we illustrate a con¬ 
nection of the conjugate FB and FF network with the directed information 
and Kramer’s causal conditioning. Let us consider the joint probability of the 
histories of x(t) and y(t), P[A(f), jV(t)]. From the definition of the conditional 
probability, we can decompose this joint probability as 


p [x(t),y(t)] = F[x(t)\y(t)]¥[y(t)} = P[y(t)\x(t)]¥[x(t)}. (37) 


However, when x{t) and y(t) are causally interacting, we can decompose the 
joint probability differently as 


F[x(t),y(t)] 


t 

n>N Xi-iM-i] 

_i —1 


t 

nw^-r] 

_i =1 


= ¥ yllx [y(t)\\x(t)}x¥ xlly [x(t)\\y(t-i)}, 


(38) 

(39) 


where P y || a ,[3^(t)||A’(t)] and ¥ x \\ y [X(t)\\y(t — 1)] are the Kramer’s causal con¬ 
ditional probability [30j. If no FB exists from y back to x, then this decom- 
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position is reduced to 


¥[x(t),y(t)] = ¥ yllx [y(t)\\x(t)} x P[ x(t)}. 


The directed information from y to x is defined as 


(40) 


I[y(t) X(t)] := ( In 




(41) 


/ p [X(t),y(t)\ 

where the joint probability of X(t) and 3^(t) is compared with the distribu¬ 


tion, P y || 2 .[jV(t)||A’(t)] x P[T(t)] 30 . Thus, I[(V(t) —> T(f)] is zero when no 


FB exists from y to x , and measures the directional flow of information from 
y back to x. In the conjugate network, the replica y' is driven only by x. 
Thus, the joint probability between X(t) and y'(t) becomes 


p [x(t),y'(t)} = p tf iix[y(t)ii*(t)]p[*(t)]. 


(42) 


Thereby, in principle^] the directed information can be calculated by obtain¬ 
ing the joint distributions, P [X{t),y{t)] and P [X(t),y'(t)], of the conjugate 
network. This relation of the conjugate network with the directed informa¬ 
tion strongly suggests that the directed information and the causal decom¬ 
position are relevant for the loop gains. Resolving this problem will lead to 
more fundamental understanding of the FB in biochemical networks because 
the directed information is found fundamental in various problems such as 


the information transmission with FB 48 , gambling with causal side infor¬ 
mation 


30 , population dynamics with environmental sensing 36 , and the 


information thermodynamics with FB 38 . This problem will be addressed 
in our future work. 
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